Quantum dynamics of an Ising spin-chain in a random transverse field 
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We consider an Ising spin-chain in a random transverse magnetic field and compute the zero 
temperature wave vector and frequency dependent dynamic structure factor numerically by using 
Jordan- Wigner transformation. Two types of distributions of magnetic fields are introduced. For 
a rectangular distribution, a dispersing branch is observed, and disorder tends to broaden the 
dispersion peak and close the excitation gap. For a binary distribution, a non-dispersing branch at 
almost zero energy is obtained. We discuss the relationship of our work to the neutron scattering 
measurement in LJH0F4. 



Calculation of real-time dynamics of a correlated quan- 
tum system with an infinite number of degrees of freedom 
are few and far between. Except for isolated examples, 
construction of real-time behavior from imaginary-time 
correlation functions (more amenable to numerical meth- 
ods) by analytic continuation is fraught with various in- 
stabilities. The theoretical challenge is particularly acute 
because neutron scattering experiments often provide a 
rather detailed map of the frequency, u), and the wave 
vector, k, dependent dynamical structure factor, 5(k, uj). 

The second motivation comes from the desire to study 
the dynamics of a quantum phase transition involving a 
zero temperature quantum critical point. In this respect, 
the Ising spin chain in a transverse fielclii2i&4 constitutes 
a schema from which much can be learned about quan- 
tum criticality^ both with and without disorder. 

The third motivation is to examine how the coherence 
of the quasiparticle excitations is modified in the presence 
of quenched disorder and is triggered by a recent neutron 
scattering experiment 6 in LiHoF4, which connects the ob- 
served low temperature behavior of S^k, uj) in terms of 
the hyperfine coupling of the electronic spins to a nu- 
clear spin bath, where the Hamiltonian of the electronic 
spins is given by an Ising model in a transverse field. Be- 
cause the hyperfine splittings are small, we could imagine 
that, on the time scale of the electronic motion, this bath 
will appear essentially quenched, modulating the quan- 
tum fluctuations characterized by the transverse field. 
This is still a bit far from the experimental system, as 
it is three-dimensional, and the Ising couplings are long- 
ranged and dipolar. Nonetheless, we shall see that there 
are tantalizing similarities between our calculated struc- 
ture factor and the experimental one in a given direction 
of the reciprocal space (2,0,0) — > (1,0,0) (in reciprocal 
units). A more appropriate comparison should be with 
quasi-one dimensional spin systems that are intentionally 
disordered. 

Finally, the role of disorder in a quantum critical sys- 
tem is an important subject in itself and is certainly not 
fully understood. In the presence of disorder, there are 
rare regions with couplings stronger than the average, 
which results in the Griffiths-McCoy singularities^. Al- 
though the effect is weak in a classical system, it be- 
comes important in quantum systems, especially in low 
dimensional^. 



The Ising model in a random transverse 
field has been studied both analytically and 
numerically^iii*i2iiiiliiSiiSiiLiSiiS and a great number 
of results of physical importance have been obtained. 
However, the dynamical structure factor S(k,uj) has 
not been computed for all k and uj in the random 
field model, although some analytical and numerical 
results are available in pure systems for special values 
of the wave vector. Here we compute the dynamic 
structure factor in the presence of two types of disorder 
distributions: a rectangular and a binary distribution. 

The one-dimensional lattice Hamiltonian we study is 



where the tr's are Pauli matrices, and J is positive and 
uniform. We shall choose the energy unit such that 
J = 1.0. The fields hi are random variables. The first 
model we study is the rectangular distribution with a 
mean h ave and a width h w . The second model is the bi- 
nary distribution in which hi is an independent random 
variable that takes two values: hs and hi, with probabil- 
ities p and (1 — p), respectively. In particular, we choose 
the parameter p to be small such that the chain is almost 
spatially homogeneous, h$ < J < hi,, to ensure that the 
system is in the paramagnetic phase. Intuitively, this 
distribution seems to capture a crude adiabatic represen- 
tation of the electronic spins coupled to a hyperfine spin 
bath discussed above, where the larger field is the applied 
transverse field. 

We first compute the time-dependent spin-spin corre- 
lation functions C(n,t) — (of (£)of +Tl ) at temperature 
T = 0, where the angular brackets denote the average 
over ground state and the overline an average over dis- 
order configurations. The algorithm for evaluating this 
quantity is similar to that described in Ref. ITR except 
that we are doing a real time calculation, which entails 
computation over complex variables. We first perform a 
Jordan- Wigner transformation^ and then cast the cor- 
relation function C{n,t) in the form of a Pfaffian, which 
is calculated efficiently, as described below; after averag- 
ing over disorder configurations, the dynamical structure 
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factor is, 



dte~ luJt C(n,t). 



(2) 



The pfaffian is the square root of the determinant of an 
antisymmetric matrix. Let X be an N x N (N is even) 
antisymmetric matrix of the form: 



X 
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where A is a 2 x 2 matrix, and £?, C are matrices of 
appropriate dimensions. From the identity: 



h 

B T A~ 1 
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where J n is the unit matrix of dimension n, we have: 

Det(X) = Det(A)Det(C + B T A~ 1 B) (5) 
Since A is also antisymmetric, of the form 



.4 : 



a 12 
-a u 



(6) 



it is easy to invert A, and hence calculate C + B T A^ 1 B. 
Since the matrix C + B T A~ X B is also antisymmetric of 
dimension (N — 2), the above procedure can be repeated, 
and the determinant of X is given by the product of 
2x2 determinants. Finally, since y/Dct(A) = a\%, the 
pfaffian of matrix X will be simply a product of N/2 
numbers obtained from those 2x2 matrices in the above 
procedure. Because it is not necessary for A to be a 2 x 2 
matrix at the upper left corner of X, for the stability of 
the algorithm, A is chosen to be a diagonal block such 
that |Det(A)| is the largest. 

The calculations were performed for a lattice of 160 
sites of which 32 sites in the middle were used for eval- 
uating C(n,t). Since we are interested in the gapped 
phase of the model where the correlation length is fi- 
nite, the size of the system was found to be sufficient by 
checking the difference between the results corresponding 
to various sizes. A window function, exp(— e\\t\ — £2^1), 
was applied to the Fourier transform in J3J) to reduce 
the cutoff effects in S(h,u))H The chosen small values 
ei = £2 = 0.05 were found to be sufficient for our nu- 
merical purposes; the range of the time integration used 
was —80 < t < 80. Free boundary condition was im- 
posed for simplicity, and the results were averaged over 
100 realizations of disorder. We explicitly checked that 
our conclusions are unaffected in going from say 20 to 
100 realizations of disorder. The reason is that we are 
interested in some robust features of the spectra of exci- 
tations in the quantum paramagnetic phase rather than 
the critical behavior or similar such subtle properties. 
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FIG. 1: (Color online) The density of states of Jordan- Wigner 
fermions, p(u)), for the rectangular distribution, for different 
disorder strengths. The p(u>) is broadened as h w increases, 
and states with gapless excitations are available at zero energy 
only in the high disorder limit. 



For Fig. 2, seeGraph2.jpg. 

FIG. 2: (Color online) The dynamic structure factor S(k,u>): 
(a)~(c) for the rectangular disorder distribution in the param- 
agnetic phase: h ave — 1.4, and h w = 0.5, 1.5,2.5 respectively; 
(d) S(k,ui) for pure system with h = 1.4 and J = 1.0. Dashed 
line shows the theoretical dispersion relation. 



To access the paramagnetic regime, in the rectangu- 
lar distribution case, not too far away from the quan- 
tum critical point, we set h ave — 1.4. The density of 
states for Jordan- Wigner fermions is shown in Fig. As 
we increase disorder, the density of states (DOS) gets 
broader, and only for strong disorder a few states ap- 
pear at zero energy, resulting in gapless excitations in 
this extreme limit. When the disorder width is small, 
for example h w = 0.5, and the system is still in the 
paramagnetic regime everywhere, there is only a single 
dispersing branch, and the excitation remains gapped, 
as in Fig. BTa). As we increase h w , the dynamic struc- 
ture factor can, in principle, have two branches because 
the DOS can acquire states at zero energy. Clearly, the 
dispersing branch, as shown in Fig. |ffib). has very low 
spectral weight everywhere. In principle, there could 
be a strong non-dispersing branch at lu = 0. However, 
since the excitation gap is closed at this value of disorder 
strength, and since both branches are significantly broad- 
ened due to randomness, the non-dispersing peak may be 
mixed with the dispersing branch, and perhaps we are 
not able to observe them separately. In the high disorder 
limit, h w = 2.5, almost all the intensity concentrates at 
(k,u>) — (0,0), but there appears be a ghost of the dis- 
persing branch, see Fig. Hfc). An important feature of 
Figs. Efb) and (c) is the horizontal stripe- like patterns, 
which indicate excitation modes that do not disperse, 
namely the localized modes. Our code was checked by 
comparing with the exact result for the pure system for 
which S(k,ui) has a single dispersive branch (see below), 



3 



though the peak is not a delta function given the finite 
size of our system; see Fig.^d). 

In contrast, the binary distribution is quite remark- 
able. We first choose hi, — 1.4 so that the system is in 
the paramagnetic regime, and set hg = 0.1, andp = 0.05. 
The density of states in Fig. shows zero energy states 
separated by a gap from the states at higher energy. 
The calculated S(k,uj) is shown in Figs. EI a) through 
(f), where hi = 1.1, 1.2, . . . , 1.6. While the dispersing 
branch is broadened due to disorder, the weight of the 
central peak around (fc, lu) = (0, 0) is so high that a non- 
dispersing branch can extend quite far away from the 
origin. 





FIG. 3: (Color online) The density of states p(w) for the 
binary distribution at hi = 1.4. Clearly, a few new states are 
allowed at zero energy, as the disorder is turned on. 



For Fig. 4, seeGraph4.jpg. 

FIG. 4: (Color online) (a)~(f) Dynamic structure factor 
S(k, uj) for the binary distribution in the paramagnetic phase. 
Parameters are J = 1.0, hs ~ 0.1, p = 0.05, and hi, = 
1.1, 1.2, . . . , 1.6 from (a) through (f). Both dispersing and 
non-dispersing branches exist, and the excitation is gapped. 

In Fig.0 we show a cut of S(k, ui) at k = 0, to illustrate 
that there is hardly any difference between the results 
averaged over 20 realizations and those averaged over 100 
realizations. 

Let us define the weight of the central peak as 



h = 



S 2 {k,uj)dkdu 



S 2 (k, uj)dkduj 



(7) 



where fl in the denominator is the entire (fc,w) domain, 
while Ail is defined as follows: we first find the max- 
imum of the central peak S' max (fc, w), and then define 
the domain Afl such that S(k,u>) > S m & x (k, lu)/\/2 for 
(k,uj) £ Aft. Note that in iQ we integrate S 2 (k,uj) 
rather than S(k,uj). This is because small numerical er- 
rors can result in slightly negative values of S(k,uj) for 
some (k,oj), which is obviously unphysical. The depen- 
dences of Ih and AO, on hi are plotted in Fig. HJ] As we 



FIG. 5: (Color online) A cut of S(0, uj) corresponding to 
Fig. Eld) for 20 and 100 realizations of disorder. 



tune hi from 1.6 to 1.1, Ih increases monotonically, while 
the region AO, shrinks. We conclude that, as the quan- 
tum phase transition is approached, the weight is trans- 
ferred from the dispersing branch to the non-dispersing 
peak. 
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FIG. 6: The binary distribution: as the quantum phase transi- 
tion is approached from the paramagnetic regime, the central 
peak intensity Ih increases monotonically, and the area Af2 
corresponding to the half the value of the peak decreases. 

In addition to the singularity due to the quantum 
phase transition, in disordered systems there is also the 
Griffiths-McCoy singularity: the disorder will drive some 
rare regions into a phase different from the rest. For the 
pure Ising chain when hi = h, the dispersion relation is£ 



2^/j 2 + h 2 -2Jhcosk. 



(8) 



The excitation gap A = 2\h — J\ occurs at k = 0. As 
h — > J, the excitation gap closes at the quantum critical 
point. However, if disorder is strong enough, there is 
nonzero probability to find some regions in which the 
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spins are strongly coupled, such that hi < J holds, and 
thus the cluster is ferromagnetic. Calculations, similar to 
that given in Ref . show that these clusters give rise 
to the non-dispersing peak at zero energy, as we describe 
below. 

For the binary distribution, the system is almost ho- 
mogeneous except for the rare regions of strongly coupled 
clusters where hi = h$ for all sites inside the cluster. At 
shorter length scales, the behavior of the pure system 
dominates, and this leads to the dispersing branch. At 
longer scales the effects of disorder become important, 
resulting in the zero energy peak. The autocorrelation 
function S(ui), which is the integral of S(k, uj) over k can 
be approximated in the following manner. The normal- 
ized probability that a given site belongs to a ferromag- 
netic cluster of length L-sites is P(L) = Lp 1 "' 1 ^ — p) 2 . 
When hi = the two-fold degenerate ground state within 
a cluster is far away from the excited states of energy of 
order 2 J; the perturbation hi = hs, will split the ground 
state by ge~ cL , where g and c are unknown positive con- 
stants determined by the details of the Hamiltonian. The 
form of the splitting results from large order in pertur- 
bation theory, however. If we treat these clusters as in- 
dependent and average over disorder, or equivalently in- 
tegrate over the cluster size L, we get 



S(uj) ~ J c\LP{L)5{u - ge~ cL ) 

(l-p) 2 hi(g/uj) g ln(p)/c 

p UJ Ul 



which diverges at ui = if 1 > |lnp|/c modulo loga- 
rithmic corrections. We have verified that indeed the 
divergence disappears for sufficiently small values of p 
(the peak of is then shifted to u> greater but close 
to zero, instead), but a detailed verification of this result 
appears to be difficult. 



It is remarkable that a numerically exact solution of 
a simplified one-dimensional model (binary distribution) 
can capture some of the experimental features in LiHoF4 
and shed interesting light on the role of disorder on the 
dynamics of a prototypical quantum critical point. It 
would be of course interesting to experimentally study 
intentionally disordered systems that are closer to the 
model discussed here. 
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